set     more off
set     matsize 800
clear   all
global  dpath = "[Userpath]\data and code"

********************************************************************************
********************************Data clean**************************************
********************************************************************************
use "$dpath\data\firmdata.dta",clear
sort countycode year

*1.County-level data
foreach var in output wastewater cod nh3n gas so2 dust {
	bysort countycode year: egen `var's = sum(`var')
	drop `var'
	rename `var's `var'
	gen l`var' = log(`var')
}
duplicates drop countycode year, force

gen llgdpC = log(lgdpC)
drop if llgdpC == .

gen llgdp = log(lgdp)
drop if llgdp == .

*2. Remove outliers
winsor2 lwastewater, replace cuts(0.5 99.5) trim
winsor2 lcod, replace cuts(0.5 99.5) trim
winsor2 lgas, replace cuts(0.5 99.5) trim
winsor2 lso2, replace cuts(0.5 99.5) trim
winsor2 ldust, replace cuts(0.5 99.5) trim

*3. Generate variables
xtset countycode year

gen manualupstream = 1 if manualyear != 33000
replace manualupstream = 0 if manualupstream == .

gen treated = 1 if autoupstream == 1 & year >= autoyear
replace treated = 0 if treated == .

gen treatedM = 1 if manualupstream == 1 & year >= manualyear
replace treatedM = 0 if treatedM == .

forvalues i=-9(1)-1 {
	local j = -`i'
	gen tN`j'= 1 if year-autoyear == `i' & autoyear != 33000
	replace tN`j' = 0 if tN`j' == .
}

forvalues i=0(1)11 {
	gen t`i' = 1 if year-autoyear == `i' & autoyear != 33000
	replace t`i' = 0 if t`i' == .
}

gen t9M  = t9+t10+t11
order t9M,before(t9) 

cd "$dpath\results"

********************************************************************************
*********************************Table******************************************
********************************************************************************

*1.Table I
reghdfe lwastewater treated , absorb(countycode citycode#year) cluster(citycode)
estimates store reg1
reghdfe lwastewater treated llgdpC, absorb(countycode citycode#year)  cluster(citycode)
estimates store reg2
reghdfe lwastewater treated treatedM llgdpC, absorb(countycode citycode#year)  cluster(citycode)
estimates store reg3
reghdfe lcod treated treatedM llgdpC, absorb(countycode citycode#year)  cluster(citycode)
estimates store reg4
esttab reg1 reg2 reg3 reg4 using TableI.rtf,replace b(%9.3f) se scalar(N r2_a) star(* 0.1 **  0.05  *** 0.01) title(Table I) keep(treated treatedM)
est clear

*2.Table II
reghdfe lgas treated treatedM llgdpC, absorb(countycode citycode#year)  cluster(citycode)
estimates store reg1
reghdfe ldust treated treatedM llgdpC, absorb(countycode citycode#year) cluster(citycode)
estimates store reg2
reghdfe lso2 treated treatedM llgdpC, absorb(countycode citycode#year) cluster(citycode)
estimates store reg3
esttab reg1 reg2 reg3 using TableII.rtf,replace b(%9.3f) se scalar(N r2_a) star(* 0.1 **  0.05  *** 0.01) title(Table II) keep(treated)
est clear

*3.Table AII
preserve
joinby citycode year using "$dpath\data\city leader age.dta",unmatched(both) 
gsort citycode secretarybirthyear
bysort citycode:g a=secretarybirthyear[_N]
drop if a==. 
drop a
g secretaryage=year-secretarybirthyear

reghdfe lwastewater treated if secretaryage<=56, absorb(countycode citycode#year) cluster(citycode)
estimates store reg1
reghdfe lwastewater treated treatedM llgdpC if secretaryage<=56, absorb(countycode citycode#year)  cluster(citycode)
estimates store reg2
reghdfe lwastewater treated if secretaryage>56, absorb(countycode citycode#year) cluster(citycode)
estimates store reg3
reghdfe lwastewater treated treatedM llgdpC if secretaryage>56, absorb(countycode citycode#year)  cluster(citycode)
estimates store reg4
esttab reg1 reg2 reg3 reg4 using TableAII.rtf,replace b(%9.3f) se scalar(N r2_a) star(* 0.1 **  0.05  *** 0.01) title(Table AII) keep(treated)
est clear
restore

*4.Table AIX Column 1
reghdfe lwastewater treated , absorb(countycode citycode#year) cluster(provincecode)
estimates store reg1
esttab reg1 using TableAIX.rtf,replace b(%9.3f) se scalar(N r2_a) star(* 0.1 **  0.05  *** 0.01) title(Table AIX Column 1) keep(treated)
est clear

*5.Table AIX Column 7
egen citycode_year=group(citycode year)
replace autoyear=0 if autoyear==33000
csdid lwastewater citycode_year,ivar(countycode) time(year) gvar(autoyear) method(drimp) cluster(provincecode) notyet agg(simple)
estimates store reg1
esttab reg1 using TableAIX.rtf,append b(%9.3f) se scalar(N r2_a) star(* 0.1 **  0.05  *** 0.01) title(Table AIX Column 7)
est clear

********************************************************************************
*********************************Figure*****************************************
********************************************************************************

*6.Figure VI
reghdfe lwastewater tN4-t8 t9M treatedM llgdp, absorb(countycode citycode#year ) 				///
        cluster(citycode)
estimates store Dynamic		
coefplot Dynamic, keep(tN4 tN3 tN2 tN1 t0 t1 t2 t3 t4 t5 t6 t7 t8 t9M) vertical ///
         recast(connect) yline(0) coeflabels(tN4 = "-4" tN3 = "-3" tN2 = "-2"   ///
		 tN1 = "-1" t0 = "0" t1 = "1"  t2 = "2" t3 = "3"  t4 = "4" t5 = "5"     ///
		 t6 = "6"  t7 = "7" t8 = "8"  t9M = "≥9")  								///
		 ytitle("Estimated Coefficients")   									///
		 xtitle("Number of years since the establishment of Automatic Stations") 			///
		 ciopts(recast(rcap) msize(medium))	
graph save "Graph" "Figure VI.gph",replace